Model Error Correction for Linear Methods of Reversible Radioligand Binding 

Measurements in PET Studies 



Hongbin Guo* ,a , Rosemary A Renaut a , Kewei Chen b , Eric M Reiman b 

a Arizona State University, Department of Mathematics and Statistics, Tempe, AZ 85287-1804- 
b Banner Alzheimer Institute and Banner Good Samaritan Positron Emission Tomography Center, Phoenix, AZ 85006 



Abstract 

Graphical analysis methods are widely used in positron emission tomography quantification because of 
their simplicity and model independence. But they may, particularly for reversible kinetics, lead to bias in 
the estimated parameters. The source of the bias is commonly attributed to noise in the data. Assuming 
a two-tissue compartmental model, we investigate the bias that originates from model error. This bias is 
an intrinsic property of the simplified linear models used for limited scan durations, and it is exaggerated 
by random noise and numerical quadrature error. Conditions are derived under which Logan's graphical 
method either over- or under-estimates the distribution volume in the noise-free case. The bias caused by 
model error is quantified analytically. The presented analysis shows that the bias of graphical methods 
is inversely proportional to the dissociation rate. Furthermore, visual examination of the linearity of the 
Logan plot is not sufficient for guaranteeing that equilibrium has been reached. A new model which retains 
the elegant properties of graphical analysis methods is presented, along with a numerical algorithm for 
its solution. We perform simulations with the fibrillar amyloid (3 radioligand [11C] benzothiazole-aniline 
using published data from the University of Pittsburgh and Rotterdam groups. The results show that 
the proposed method significantly reduces the bias due to model error. Moreover, the results for data 
acquired over a 70 minutes scan duration are at least as good as those obtained using existing methods 
for data acquired over a 90 minutes scan duration. 

Key words: Bias; graphical analysis; Logan plot; PET quantification; PIB; Alzheimer's disease; 

distribution volume. 
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1. Introduction 



Graphical analysis (GA) has been routinely used for quantification of positron emission tomography 
(PET) radioligand measurements. The first GA method for measuring primarily tracer uptakes for 
irreversible kinetics was introduced by Patlak, G], 0]> an d extended for measuring tracer distribution 
(accumulation) in reversible systems by Logan, [3|. These techniques have been utilized both with input 
data acquired from plasma measurements and using the time activity curve from a reference brain region. 
They have been used for calculation of tracer uptake rates, absolute distribution volumes (DV) and DV 
ratios (DVR), or, equivalently, for absolute and relative binding potentials (BP). They are widely used 
because of their inherent simplicity and general applicability regardless of the specific compartmental 
model. 
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The well-known bias, particularly for reversible kinetics, in parameters estimated by GA is commonly 
attributed to noise in the data, [1, IB, 0] , and therefore techniques to reduce the bias have concentrated 
on limiting the impact of the noise. These include (i) rearrangement of the underlying system of linear 
equations so as to reduce the impact of noise yielding the so-called multi-linear method (MAI), [H], and 
a second multi-linear approach (MA2) , 0j , (h) preprocessingusing the method of generalized linear least 
squares (GLLS ), [a] , yielding a hybrid GLLS-GA method, [21], fiii) use of the method of perpendicular 
least squares, [101 ]. also known as total least squares (TLS), 11 111, (iv) likelihood estimation, 12|, (v) 



Tikhonov regularization [13], (vi) principal component analysis, [141 ]. and (vii) reformulating the method 



of Logan so as to reduce the noise in the denominator, 15[]. Here, we turn our attention to another 
important source of the bias: the model error which is implicit in GA approaches. 

The bias associated with GA approaches has, we believe, three possible sources. The bias arising due to 
random noise is most often discussed, but errors may also be attributed to the use of numerical quadrature 
and an approximation of the underlying compartmental model. It is demonstrated in Section [2] that not 
only is bias an intrinsic property of the linear model for limited scan durations, which is exaggerated by 
noise, but also that it may be dominated by the effects of the model error. Indeed, numerical simulations, 
presented in Section [U demonstrate that large bias can result even in the noise-free case. Conditions 
for over- or under-estimation of the DV due to model error and the extent of bias of the Logan plot 
are quantified analytically. These lead to the design of a bias correction method, Section O which still 
maintains the elegant simplicity of GA approaches. This bias reduction is achieved by the introduction of 
a simple nonlinear term in the model. While this approach adds some moderate computational expense, 
simulations reported in Section I4.3| f or the fibrillar amyloid (3 radioligand [11C] benzothiazole-aniline 
(Pittsburgh Compound- B [PIB]), [la ], illustrate that it greatly reduces bias. Relevant observations are 
discussed in Section [5] and conclusions presented in Section [6l 

2. Theory 

2.1. Existing linear methods 

For the measurement of DV, existing linear quantification methods for reversible radiotracers with a 
known input function, i.e. the unmetabolized tracer concentration in plasma, are based on the following 
linear approximation of the true kinetics developed by Logan, [jj: 

MAO: I C T (r)dr«DV f C p (r)dr - bC T {t). (1) 
Jo Jo 

Here Cr(i) is the measured tissue time activity curve (TTAC), C p (t) is the input function, DV represents 
the distribution volume and quantity b is a constant. With known Ct(£) and C p (t) we can solve for DV 
and b by the method of linear least squares. This model, which we denote by MAO to distinguish it 
from MAI and MA2 introduced in [H], approximately describes tracer behavior at equilibrium. Dividing 
through by Cr(t), showing that the DV is the linear slope and —b the intercept, yields the original Logan 
graphical analysis model, denoted here by Logan-GA, 

Logan - GA : K Dv _ 6> (2) 

in which the DV and intercept —b are obtained by using linear least squares (LS) for the sampled version 
of ([2]). Although it is well-known that this model often leads to under-estimation of the DV it is still 
widely used in PET studies. An alternative formulation based on ([I]) is the so-called MAI, 

DV [ l 1 f* 

MAI : C T (i) **— J C p (r)dr --J C T (r)dr, (3) 



for which the DV can again be obtained using LS [a*] ■ Recently another formulation, obtained by division 
in ([1]) by C p (t) instead of Cr(i), has been developed by Zhou et al, [la]. But, as noted by Varga et al in 
[lo| the noise appears in both the independent and dependent variables in ((2J) and thus TLS may be a 
more appropriate model than LS for obtaining the DV. Whereas it has been concluded through numerical 
experiments for tracer [ 18 F]FCWAY and [ n C]MDL 100,907, 5|, that MAI © performs better than other 
linear methods, including Logan-GA ([2J, TLS and MA2 [7J, la], none of these techniques explicitly deals 
with the inherent error due to the assumption of model MAO (jlj) . The focus here is thus examination of 
the model error specifically for Logan-GA and MAI, from which a new method for reduction of model 
error is designed. 



2.2. Model error analysis 

The general three-tissue compartmental model for the reversible radioligand binding kinetics of a 
given brain region or a voxel can be illustrated as follows, 17J, 1 181 ] : 



C p (t) 



C F (t) 



C NS it) 



Cs(t) 



Figure 1: Three-tissue compartmental model of reversible radioligand binding dynamics. 



Here C p (t) (kBq/ml) is the input function, i.e. the unmetabolized radiotracer concentration in plasma, 
and Cp(t), Cns(£) and Cs(t) (kBq/g) are free radioactivity, nonspecific bound and specific bound tracer 
concentrations, resp., and K\ (ml/min/g) and ki (1/min), % = 2, • • • ,6, are rate constants. The DV is 



related to the rate constants as follows [l9J], 



Ki k-t kx. . 

DV = ^ (1 + H + ^»' (4) 

The numerical implementation for estimating the unknown rate constants of the differential system 



illustrated in Figure[T]is difficult because three exponentials are involved in the solution of this system, [181 ] . 
Specifically, without the inclusion of additional prior knowledge, the rate constants may be unidentifiable, 
[201 ] . Fortunately, for most tracers it can safely be assumed that Cns an d Cf reach equilibrium rapidly 
for specific binding regions. Then it is appropriate to use a two-tissue four-parameter (2T-4k) model 
by binning Cns(*) and Cp(t) to one compartment Cp+Ns(t) = Cp(t) + Cns(*)- This is equivalent to 
taking k§ = k$ = 0, and hence Cns(*) = 0. On the other hand, for regions without specific binding 
activity, we know C$(t) = which is equivalent to taking k% = k^ = 0, and it is again appropriate for 
most radioligands to bin Cns(^) an d Cp(t). The one-tissue compartmental model is then appropriate for 
regions without specific binding activity. For some tracers, however, for example the modeling of PIB in 
the cerebellar reference region, the best data fitting is obtained by using the 2T-4k model without binning 
Cns (t) and Cp (t) , (2l| . Assuming the latter , the DV is given by K\ jk^ ( 1 + ^3 /k± ) , and K\ / ki ( 1 + k^ / k§ ) , 
for regions with and without specific binding activity, resp. Ignoring the notational differences between 
the two models, for regions with and without specific binding activity, they are both described by the 
same abstract mathematical 2T-4k model equations. Here, without loss of generality, we present the 
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2T-4k model equations for specific binding regions, 



' t} = ^C^-ih + k^Cp+Ns^ + hCsit) (5) 



dt 
dCg(t) 

dt 



k 3 C F+NS (t) - hC s (t). (6) 



To obtain the equations appropriate for regions without specific binding activity, Cs(t) is replaced by 
Cns(*) and &3 and k 4 are interpreted as the association and dissociation parameters of regions without 
specific binding activity. To simplify the explanation C$(t), k$ and k 4 are used throughout for both 
regions with and without specific binding activity, with the assumption that Cs(t), k% and k 4 should 
automatically be replaced by Cns(*) 5 &5 and k§ respectively, when relevant. 
The solution of the linear differential system dS])-© is given by 

C F+NS (t) = (aie- ait + 6 ie - Q2t )®C P (t) (7) 
C s (t) = a 2 (e- Ql '-e- a2 ')®C p (t) (8) 

where <g) represents the convolution operation, 

"1,2 = (k 2 + k 3 + k 4 =F + k 3 + k 4 ) 2 - Ak 2 k 4 )/2, and 

Ki(/c 4 -ai) Ki(a 2 -k 4 ) K x k z 

a\ = , b\ = , &nda 2 = . (9) 

a 2 — a\ ol 2 — ol\ a 2 — ai 

The overall concentration of radioactivity is 

C T (t) = C F+NS (t) + C s (t) = ((ai + a 2 ) e ~ Qli + (6i - a 2 )e~ a ^) ® C p (t). (10) 
Integrating and rearranging yields 

/W)dr = DV f C p (r)dr - ^±C F+NS (t) - ** \ k \ + h C s (t), (11) 
Jo Jo kzki k 2 k 4 

= DV^ CpW d T -^lc T ( t )-ic s(t ). (12) 

This is model ([T]) when Cs(t) is linearly proportional to Ct(^) for a time window within the scan duration 
of T minutes. The accuracy of linear methods based on ([1]) is thus dependent on the validity of the 
assumption that Cs(t) and C F +Ns(t) are approximately linearly proportional to Cr(t) over a time window 
within [0, T]. Logan observed that C F+ Ns(t) and Cs(t) are roughly proportional to Cr(i), after some 
time point t* , [1] . If the assumption of linear proportionality breaks down for the given window, [t* , T] , 
bias in the estimated uptake rate or DV will be introduced, as shown later in Section 14.31 due to the 
intrinsic model error of a GA method. Indeed, in Section 15.11 we show that, for the PIB radioligand 
on some regions with small k 4 , there is no window within a 90 minutes scan duration where Cs(t) and 
Cr(i) are linearly proportional. This is despite the apparent good linearity, visually, of the Logan plot 
of J Ct?{t)&t / C^it) against J* C p (r)dr/CT(i). Waiting for equilibrium, which may take several hours, 
is impractical in terms of patient comfort, cost and measurement of radioactivities. 

The limitation of the constant approximation can be analysed theoretically. Because a 2 >> ot\ > 
and C p (t) is very small for large time the convolution e~ a2t <8> C p (t) = J * e~ a2 ^~ r ^Cp(r)dT is relatively 
small. We can safely assume that the ratio of e~ a2t <8> C p (t) to e~ Ql * <g> C p (t) is roughly for t > t*. Then 
Cs(t), see equation (jSJ), is approximately proportional to e~ ait (8) C p (t) for t > t* . In our tests with PIB, 
the neglected component a 2 e~ a2t <8> C p (t) is less than 8%Cs(t) for t > 35 min.. On the other hand, this 
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is not the case for Cp+Ns{t), see equation ((7J), because a\ and b\ need not be of the same scale. For 
example, if fc 4 « k 2 + £3 we know b\/a\ ~ (/c 2 + ^3)/ (2fc 4 ) from ([9]), thus 61 >> ai > 0. Specifically, 
b\e~ a2t ®C p (i) may not be small in relation to a\e~ ait <8>C p (t). Thus, it is not appropriate, as is assumed 
for the Logan-GA ([2|) and other linear methods derived from MAO, to approximate 

-( f \ _ h + h C F+NS (t) fc 2 + fc 3 + fc 4 Cs(t) , , 

S[) ~ k 2 k A ' C T (t) + fc 2 fc 4 'C T (i)' 1 j 

as constant for t G One may argue that if (ai + a 2 )/(bi — a 2 ) is close to 1 the term e~ a2t ® C p (t) 

in Cr(t) could be ignored. Then the ratio of Cr(t) to Cs(t) would be close to constant after t*, and the 
resulting estimates of the DV using Logan-GA ([5]) and MAI ([3]) would be reasonable. While it is easy 
to verify that (a± + a 2 )/(6i — a 2 ) is positive and bounded above by one, this fraction need not be close 
to its upper bound. Indeed, for realistic test data, see Table [JJ 0.05 < (ai + a 2 )/(6i — a 2 ) < 0.65. The 
simulations presented in Tables [2] and [3] validate that a small value of this fraction may cause a problem 
in the estimation of the DV using the linear Logan-GA and MAI methods. 

It is immediate using Cx(i) = Cp+Ns(t) + Cs(t), and positivity of both CF+Ns(t)/CT(t) and 
Cs(t)/Cr(t), that s(t) is bounded above and below, 

kz + k A ^ ^ k 2 + k 3 + fc 4 _ k 3 + fc 4 + J_ 

/c 2 &4 fc 2 fc 4 fc 2 /c 4 fc 4 ' 

and 1/&4 determines the variation in s(t). If /c 4 is small the bound is not tight and the DV estimated by 
Logan-GA, or a linear method derived from MAO, may not be accurate, see for example the regions of 
interest (ROIs) 1, 3 and 6 in the test examples reported in Table [TJ We reiterate that, by the discussion 
above, the variation for ROI 6, within which no specific binding activity exist, is determined by 1/&6- 
This relationship between the size of k^ and the bias in the Logan-GA estimate of the DV is illustrated 
in Figure [8] of Section 15.11 for the test data of Table [TJ 

2. 3. Model error of Logan equation 

The complete mathematical result for the model error of Logan-GA and MAO is presented in the 
Appendix. Similar results, omitted here to save space, can be obtained for MAI. The main conclusion 
is that both Logan-GA and MAO can lead to an over-estimation of the DV. This contrasts the standard 
view of these methods. We summarize in the following theorem, for which the main idea is to show that 
replacing (|13|) which occurs on the right hand side of (jlip by a constant intercept b introduces an error 
in the least squares solution for the DV which can be specifically quantified. 

Corollary 1. Suppose Logan-GA, or respectively MAO, are used for noise- free data acquired for n frames 
with frame time ti, i = 1, • • ■ , n and t* = t[. Then, with s(t) as defined in \13\) , for each method the same 
conclusions are reached: 

• The DV is over-estimated (under- estimated) if s(t),t G [tj,£ra] 7 a non-constant decreasing (in- 
creasing) function, and 

• the DV is exact ifs(t),t G is o> constant function; 

Let DVt be the true value of the DV, and define the variation of a function over [tj,tn] by 

V(x(t)) = \ max x(i) - min x(t)|. (15) 

te[t u t n ] te[ti,tn] 
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Then the bias in DVl calculated by Logan- G A is bounded by 



(n- J + 1)^ pi 



|DV L - DV T | < V{a(t)), (16) 

2_, (Pi - Pj) 

w/iere p, = J^' Cp(r)dT/Cr(*i)- 

This theorem is an immediate result of Lemma [3] and Corollary [3] in the Appendix for the vectors 
obtained from the sampling of the functions 

s(t) = C F+NS (t) + — C s (t), and 

r(*) = fcri^dr, p(t) = f C p (t)&t, q(t) = C T (i), 
J •/ 

at discrete time points t = ti, ■ ■ ■ ,t n . The relevant vectors are defined by r = r/q, p = p/q, s = s/q, 
where the division corresponds to componentwise division. It is easy to check that all these vectors are 
positive vectors, p, p, r and r are non-constant increasing vectors and q is decreasing. Thus all conditions 
for Lemma [3] and Corollary [3] are satisfied. Note that in the denominator of (]16p the simplification 

(n — l + 1) J27=i(Pi) 2 ~ (Yn=i Pi) 2 = ^i^j.i<i,j<n(Pi ~ Pj) 2 is nse( ^- In tne l atter discussion we may use 
the variation (increasing or decreasing) of Cs(t)/Cr(t) instead of that of s(t) because 

It is not surprising that the properties of Logan-GA and MAO are similar. Indeed, MAO is none 
other than weighted Logan-GA with weights Cr(tj), which changes the noise structure in the variables. 
In contrast to the conventional under-estimation observations, it is suprising that the DV may be over- 
estimated. However, the over-estimation is indeed observed in the tests presented in Section [4.21 and [4.31 
Inequality (I16p indicates that Logan-type linear methods will work well for data for which V(s) is flat. 
Unfortunately, V(s) may become flat only for a late time interval. Thus our interest, in Section [3l is to 
better estimate the DV using a reasonable (practical) time window, which may include the window over 
which Cs(t)/Cr(i) is still increasing. Our initial focus is on the modification of Logan-type methods. 
Then, in Section [J] we present numerical simulations using noise-free data which illustrate the difficulties 
with Logan-GA and MAI, and support the results of Theorem [H 



3. Methods 

In the previous discussion we have seen the theoretical limitations of the Logan-GA and MAI methods. 
Here we present a new model and associated algorithm which assists with reducing the bias in the 
estimation of the DV. 

Observe that, «2 >> cei, implies that Cs = a2e~ ait <S> C p (t) + e(i), where e(t) can be ignored for 
t > t* . Therefore, for t > t* (|12p can be approximated by a new model as follows 

t ft 
C T (T)dT^DV C p (T)dT-AC T (t)-Be- ait ®C p (t), (17) 
Jo 

where A = (ks + k/^jh^k.^ and B = aijk^. This suggests new algorithms should be developed for 
estimation of parameters DV, A, B and a\. Here, a new approach, based on the basis function method 
(BFM) in [22|, in which oc\ is discretized, is given by the following Algorithm. 
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Algorithm 1. Given C p (ti) and Ct(^) for i = 1, • • • , n and t* = ti, the DV is estimated by performing 
the following steps. 

1. Calculate DV and intercept —b, using Logan-GA. 

2. Set af in = 0.001 and af ax = min(l, 2/6) if b > otherwise af ax = 1. 

(7) 

3. Form discretization , j = 1 : 100 for ct\, with equal spacing logarithmically between af m and 

„,max 

4. For each j solve the linear LS problem, i.e. cast it as a multiple linear regression problem with 
Jq Ct{t) dr as the dependent variable. 

DV [ C p (r)dT - AC T (t) - B I e~ a ^ T C p (t - r)dr « I C T (r)dr (18) 
Jo Jo Jo 

with data at ti, i = I, - ■ ■ ,n, to give values BY^, A® and B^.: 

5. Determine af ^ for which residual is minimum over all j. Set DV, A and B to be DV *), A^ 
and B^ ' , resp. 

Remarks: 



The interval for a\ is determined as follows: First the lower bound 0.001 for a\ is suitable for most 
tracers, but could be reduced appropriately. This lower bound is not the same as that on used 



in BFM, in which 9 is required to be greater than the decay constant of the isotope, 22|. Second 
by point ([2]) of Corollary [3] in the Appendix A, b should be positive and near the average value of 
s(t), where, by (H}, ^ < §(t) < ^gf±^ ■ On the other hand, k2+ k %+ kA « ^ if ±k 2 k 4 is small 
relative to (k 2 + ^3 + k^) 2 . Thus, ct\ is linked with b through s(t). This is used to give the estimate 
of the upper bound on a±. Practically, it is possible that the Logan-GA may yield an intercept 
b < 0, then we set a™ ax = 1. 

2. Numerically, because Jq C p (r)dr is much larger than both Cr(i) and Cs(t) for t > t*, the estimate 
of DV is much more robust to noise in the formulation, including both model and random noise 
effects, than are the estimates of A and B. Therefore, while A and B may not be good estimates 
of {kz + k\)j{k,2k\) and a2/^4, resp. for noisy data, the estimate of DV will still be acceptable. 
Consequently, it is possible that Logan-GA and MAO will produce reasonable estimates for DV, 
even when the model error is non negligible. 

3. The algorithm can be accelerated by employing a coarse-to-fine multigrid strategy. The coarser 
level grid provides bounds for the fine level grid. The grid resolution can be gradually refined until 
the required accuracy is satisfied. 

4. Experimental Results 

We present a series of simulations which first validate the theoretical analysis of Section [2] for noise- free 
data, and then numerical experiments which contrast the performance of Algorithm [1] with Logan-GA, 
MAI and nonlinear kinetic analysis (KA) algorithms for noisy data. 

4-1. Simulated Noise-Free Data 

We assume the radioligand binding system is well modeled by the 2T-4k compartmental model and 
focus the analysis on the bias in the estimated DV which can be attributed to the simplification of the 
2T-4k model. For the simulation we use representative kinetic parameters for brain studies with the PIB 



tracer. These kinetic parameters, detailed in Tabled! are adopted from published clinical data, [2ll. 1 2 
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The simulated regions include the posterior cingulate (PCG), cerebellum (Cere) and a combination of 
cortical regions (Cort). The kinetic parameters of each ROI are also associated with the subject medical 
condition, namely normal controls (NC) and Alzheimer's Disease (AD) diagnosed subjects. The kinetic 
parameters for the first seven ROIs are from [21[ while the last four are from 23f] . Rate constants for 
ROIs 5 to 11 are directly adopted from the published literature, while those for ROIs 1 to 4 are rebuilt 
from information provided in 21[. The values for ROIs 1 to 4 and 8 to 11 represent average values for 
each group, while those for ROIs 5 and 6 are derived from one AD subject and those for ROI 7 from 
another AD subject. 



Table 1: Rate constants for eleven ROIs, including PCG, Cere, and Cort, for AD and NC adopted from [ill. |23|. For ROIs 6, 
7, 10 and 11 no specific binding activity is assumed, i.e. k$ = k± — 0, DV = Ki / k2(l + k$ / k§); while for ROIs 1 to 5, 8 and 9 
we assume that the free and nonspecific compartments rapidly reach equilibrium, i.e. k$ = k& = 0, DV = K±/k2(l + kz/ki). 
Coefficients ai, foi and 02 are defined in ||9]|. The values for ROIs 1 to 4 and 8 to 11 represent average values for each group, 
while those for ROIs 5 and 6 are derived from one AD subject and those for ROI 7 from another AD subject. 



ROI/Group 


Area 


Kt 


k 2 


^3 




h 




DV 


0,1+0,2 

61 -a 2 


1/NC 


Cort 


0.250 


0.152 


0.015 


0.0106 








3.9722 


0.11 


2/ AD 


Cort 


0.220 


0.113 


0.056 


0.023 








6.6872 


0.65 


3/NC 


PCG 


0.250 


0.150 


0.015 


0.0106 








4.0252 


0.11 


4/ AD 


PCG 


0.220 


0.100 


0.050 


0.017 








8.6706 


0.63 


5/ AD 


PCG 


0.262 


0.121 


0.044 


0.015 








8.5168 


0.44 


6/ AD 


Cere 


0.273 


0.144 








0.007 


0.005 


4.5500 


0.05 


7/AD 


Cere 


0.333 


0.172 








0.029 


0.042 


3.2728 


0.26 


8/NC 


Cort 


0.250 


0.140 


0.020 


0.018 








3.7480 


0.18 


9/AD 


Cort 


0.220 


0.110 


0.050 


0.025 








5.9841 


0.63 


10/NC 


Cere 


0.270 


0.140 








0.020 


0.026 


3.4353 


0.20 


11/ AD 


Cere 


0.260 


0.130 








0.020 


0.025 


3.5810 


0.22 



The noise-free decay-corrected input function is adapted from the plasma measurements for a NC 
subject as presented in Figure 3(A) of 2l|]. Using the data from that figure we convert to kBq/ml 
under the assumption of a 100kg body mass, and obtain the functional representation for C p (t) = u(t), 
(kBq/ml), which is illustrated in Figure [2j 



u(t) 





407.4933(4 - 0.3) 
-436.64 + 384.208 

46.6747(4 + 0.24)" 2 - 2560 + 5.7173(4 + 0.24)" - 5644 



* € [0,0.3] 
4 € [0.3,0.6] 
4 € [0.6,0.76] 
4 > 0.76. 



(19) 



Using this input function and the eleven data sets given in Table Q] eleven noise-free TTACS, Ct(4) 
(kBq/ml), are generated using the 2T-4k model. The scanning protocol, consistent with that adopted in 
[211 ]. has frame durations, A4j, measured in minutes, 4 x 0.25, 8 x 0.5, 9x1, 2x3, 8x5 and 3 x 10. 
The last eight frames, which fall in the window from 35 to 90 minutes, are chosen for the time window 
over which we assume that equilibrium is achieved. A scan duration of 90 minutes is common for most 
PIB-PET dynamic studies, 



4-2. Examples of over- estimation for Logan- G A and MAI 

Theorem [T] predicts that the DV will be over-estimated when s decreases. This is validated for data for 
the simulated ROIs. The estimates of the DV, for scan durations T = 90 minutes with t* = 35 minutes, 
and T = 240 minutes with t* = 100 minutes, are reported in Table [2j The extended time window 
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♦ with noise 

°5 10 20 30 40 50 60 70 80 90 

time (min) 

Figure 2: The true input function as given by (|19p . and the simulated measurements with noise. The simulated measurements 
are generated by (|21[) with CVs = 0.05, e = 50%, /j, = 0.5ml and Awi — 100 seconds. The function over the initial 5 minutes 
is illustrated in the inset. 



is generated by adding 15 frames each of 10 minutes length. Indeed, the over-estimation predicted in 
Theorem [1] is confirmed for ROI 7, for which the decrease of Cs(£)/Cr(£) and, hence s after 35 minutes, 
is clearly illustrated in Figure [6J Moreover, Cs(£)/Cr(i) is decreasing after 100 minutes for all ROIs 
except ROI 6, see Figure |6(b)| and in all but this case the values of DV are over-estimated. We note 
that s is nearly flat on the selected windows, [i*,T] for the cases in which the over-estimation of DV is 
small. These results further validate the conclusions of Theorem [TJ Additionally, the use of the long scan 
duration of 240 minutes leads to estimates with less overall bias because the variation in Cs(t)/CT(t) 
is smaller over [100mm., 240mm.] than over the earlier window. Equivalently, as given by (I16p . a small 
variation in s guarantees a small error in the estimated DV. Clearly, linear methods based on the MAO 
model work well during the equilibrium phase. Unfortunately, this equilibrium may be reached too late 
for practical application, see for example ROI 6 in Figure [o(b)[ for which approximate equilibrium is not 
reached until 3 hours. The results with 90 minutes scan duration show that better estimates are obtained 
for larger (a± + 02) /{b\ — 02), which consistently supports the analysis in Section [231 

In these simulations the accurate data and integrals are used so as to assure that the results are not 
impacted by use of a low accuracy numerical quadrature but instead are focused on the effects of the 
model error of Logan-GA and MAI. It is interesting to note, however, that the error introduced by the 
numerical quadrature always lowers the estimate of the DV, see Section 15.21 Moreover, the noise from 
other sources may have a similar impact. This is a topic for future research. 

4-3. Algorithm Performance for Noise-Free Data 

We contrast the performance of Algorithm [T] with Logan-GA, MAI and KA for noise-free data. The 
use of a long scan duration (up to 90 minutes) is to assure that equilibrium is achieved as needed for GA 
methods. For a method for which the bias due to model error is not impacted by the need for equilibrium, 
a shorter scan duration is preferred. For the results presented in Table [3] the DV is calculated for the 
noise-free case over a scan duration of just 70 minutes with t* = 35 minutes. Accurate integrals are used 
so as to focus the conclusions on the impact of the model error. 

The KA solutions were obtained using two different optimization algorithms for the solution of the 
highly nonlinear problem, the interior point and the Marquardt-Levenberg methods, Matlab@ functions 
fmincon and lsqnonlin, resp. In order to provide the most fair comparison the results presented are 
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Table 2: The DV calculated using Logan-GA and MAI with noise-free data and accurate integrals. DV is calculated for 
scan durations T = 90 minutes with t* = 35 minutes, and T = 240 minutes with t* = 100 minutes. The percentage bias is 
listed in parentheses. 



ROI 
ID 


True 
DV 


35-90 min 


100-240 min 


Logan-GA 


MAI 


Logan-GA 


MAI 


1 


3.9722 


3.549(- 10.65%) 


3.542(-10.84%) 


3.981(0.22%) 


3.977(0.12%) 


2 


6.6872 


6.585(-1.53%) 


6.577(-1.65%) 


6.709(0.33%) 


6.709(0.33%) 


3 


4.0252 


3.593(-10.73%) 


3.586(-10.92%) 


4.034(0.22%) 


4.030(0.11%) 


4 


8.6706 


8.342(-3.79%) 


8.331(-3.92%) 


8.687(0.19%) 


8.685(0.16%) 


5 


8.5168 


8.129(-4.55%) 


8.117(-4.69%) 


8.536(0.23%) 


8.533(0.19%) 


6 


4.5500 


3.204(-29.58%) 


3.208(-29.50%) 


4.281(-5.91%) 


4.273(-6.10%) 


7 


3.2728 


3.300(0.82%) 


3.298(0.76%) 


3.286(0.41%) 


3.288(0.45%) 


8 


3.7480 


3.635(-3.01%) 


3.625(-3.28%) 


3.780(0.84%) 


3.779(0.84%) 


9 


5.9841 


5.910(-1.23%) 


5.902(-1.37%) 


6.007(0.38%) 


6.007(0.39%) 


10 


3.4353 


3.416(-0.57%) 


3.408(-0.78%) 


3.462(0.77%) 


3.463(0.80%) 


11 


3.5810 


3.552(-0.81%) 


3.544(-1.04%) 


3.608(0.75%) 


3.609(0.79%) 



for fmincon, which gave the better solutions. The KA solution is very dependent on provision of a good 
initial value. If the initial values of k% and k^ are taken very close to their true values, the estimate of the 
DV may be nearly perfect. Here we use initial values for K±, &2, k% and k^ set to [0.2, 0.1, 0.01, 0.001]. 

For Logan-GA and MAI, solutions were also calculated for the scan duration of T = 90 minutes with 
t* = 35 minutes as illustrated in Table [5J The KA results, not given, which do not require the attainment 
of equilibrium were comparable for both scan durations as expected. This independence with respect to 
the requirement of attainment of equilibrium was also observed for Algorithm [T] except for ROI 6. In this 
case the neglected part in model (|L7j) is relatively large as compared to that for the other ROIs, i.e. the 
ratio of e~ a2t <X> C p (t) to e~ ait <g) C p (t) for ROI 6 is greater than that for the other ROIs. A significant 
reduction in the bias for ROI 6 from —12.71% (70 min.) to —7.39% (90 min.) was observed. It is clear, 
by comparing the results with those in Table [21 that Algorithm [1] for a scan duration of just 70 minutes 
is much more accurate for the calculation of the DV than are Logan-GA and MAI using scan durations 
of 90 minutes. 

Table 3: DV calculated by Logan-GA, MAI, KA and Algorithm Q] for a 70 minutes scan duration with t* = 35 minutes. In 
each case the percentage bias is listed in parentheses. 



ROI 


Logan-GA 


MAI 


KA 


Algorithm [T] 


1 


3.395(- 14.54%) 


3.392(-14.61%) 


3.928(-1.12%) 


4.014(1.05%) 


2 


6.511(-2.64%) 


6.506(-2.71%) 


6.552(-2.02%) 


6.777(1.34%) 


3 


3.436(-14.65%) 


3.433(-14.71%) 


3.982(-1.08%) 


4.066(1.00%) 


4 


8.163(-5.86%) 


8.157(-5.92%) 


8.535(-1.56%) 


8.743(0.83%) 


5 


7.931(-6.88%) 


7.925(-6.95%) 


8.383(-1.57%) 


8.530(0.16%) 


6 


3.004(-33.97%) 


3.007(-33.92%) 


4.675(2.74%) 


3.972(-12.71%) 


7 


3.293(0.63%) 


3.292(0.58%) 


3.188(-2.59%) 


3.277(0.12%) 


8 


3.555(-5.15%) 


3.549(-5.30%) 


3.679(-1.84%) 


3.784(0.95%) 


9 


5.847(-2.28%) 


5.842(-2.37%) 


5.859(-2.10%) 


6.008(0.40%) 


10 


3.376(-1.73%) 


3.371(-1.87%) 


3.361(-2.17%) 


3.451(0.47%) 


11 


3.506(-2.09%) 


3.501(-2.24%) 


3.505(-2.11%) 


3.585(0.10%) 



In contrasting the results with respect to only the bias in the calculation of the DV it is clear that 
Algorithm [T] leads to significantly more robust solutions than Logan-GAl and MAI for noise-free data. 



10 




10" 1 10 1 1 10 2 20 4 6 80 1 00 



time (min) time (min) 

(a) (b) 

Figure 3: (a) The coefficients of variation CVt for the noisy TTAC associated with ROf 3, obtained with Sc = 1 and Sc — 2, 
resp. and CVr, for the input function calculated for e = 50%, /j, = 0.5ml and At«i = 100 seconds, (b) The noise-free and 
noisy TTACs for ROI 3 obtained with Sc — 1 and Sc — 2, resp. 



On the other hand, the KA approach can lead to very good solutions, comparable and perhaps marginally 
better than Algorithm [TJ For ROI 6, for which the KA solution is significantly better, we recall that the 
solution depends on the initial values of the parameters. Changing the initial k§ to 0.01, the resulting 
bias in the DV of ROI 6 calculated by KA is increased to 31.75%. On the other hand, Algorithm [TJ is 
not dependent on specifying initial values, and is thus more computationally robust. 

4-4- Experimental Design for Noisy Data 

While the results with noise-free data support the use of Algorithm [lj it is more critical to assess its 
performance for noise-contaminated simulations. The experimental evaluation for noisy data is based on 
the noise-free input u(t) and noise-free output Cr(t), one output TTAC for each of the eleven parameter 
sets given in Table [TJ Noise contamination of the input function and these TTACs is obtained as follows. 

4-4-1- The Noise- Contaminated TTAC Data 

For a given noise-free decay-corrected concentration TTAC, Cr(i), Gaussian (G(0, cr(CT(t))) noise at 
each time point t; t is modeled using the approach in 0, Eol . 0). The standard deviation in the noise at 
each time point tj, depends on the frame time interval Ati in seconds, the tracer decay constant A (0.034 
for 11 C) and a scale factor Sc 

a(C T (U)) = sJ CTi ^ Mi . (20) 

The resulting coefficients of variation CVt (ratio cr(Cr(ii)) to Cr(ti))' for scale factors 1 and 2, are 
illustrated in Figure [3l 

4- 4-2- The Noise-Contaminated Input Function 

The noise in the input function can be attributed to two sources, system and random noise. Although 
the random 7-ray emission follows a Poisson distribution, we use the limiting result that a large mean 
Poisson distribution is approximately Gaussian to model this randomness as Gaussian. Thus both sources 
are modeled as Gaussian but with different variance. Consider first the following model for determining 
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the randomness of the 7— ray emissions. Suppose a fi ml blood sample is placed in a 7-ray well counter 
which has efficiency e and the measured counts over Awi seconds are n(t,). Then the measured decay 
corrected concentration (kBq/ml) is 

n{tj)e xt i 
Up[ti) ~ lOOOA^/xe' 

where 1000 is a normalization factor to convert the counts to "kilo" counts. Then, assuming that the 
mean of C p (ti) (or its true value) is u{tj) as given in (|19p . the standard deviation in the measurement 
of Cp(ij) due to random effects is <7r(C p (^)) = 1/ u(ti)e xt * / (lOOOAuij^e). The coefficient of variation, 
CVr = 0R(C p (ij)) /u(ti) , which results from this random noise is shown in Figure 02 It is assumed in the 
experiments that each blood sample has volume fi = 0.5ml, the count duration is Aw^ = 100 seconds and 
the well counter efficiency is e = 50%. Then, denoting the coefficient of variation due to system noise by 
CVs, the noise-contaminated input is given by 

C V {U) = u{U){\ + (CV R + CVs)??,), (21) 

where m is selected from a standard normal distribution (G(0, 1)), and in the simulations we use CVs = 
0.05, see Figure El 



4-5. Experimental Results for Noisy Data 

Two hundred random noise realizations are generated for each input-TTAC pair, and for each noise 
level (Sc = 1, 2). The distribution volume is calculated for each experimental pair using Logan-GA, MAI, 
KA and Algorithm [TJ In each case two scan durations are considered, 70 and 90 minutes respectively, 
and t* = 35 minutes. Unlike the noise-free case, the numerical quadrature for J *C p (r)dr uses only the 
samples at scan points C p (ij). 

We present histograms for the percentage relative error of the bias 100(DV es t — DVt)/DVt in order 
to provide a comprehensive contrast of the methods. Figure H] shows the histograms for all eleven ROIs, 
with the range of the error for each method indicated in the legend. The figures (a)-(b) are for scan 
windows of 90 minutes, for noise scale factors Sc = 1 and Sc = 2 while (c)-(d) are for scan windows of 70 
minutes. Figure [5] provides equivalent information for a representative cortical region ROI 3. It is clear 
that the distributions of the relative errors for KA and MAI are far from normal; KA has a significant 
positive tail while Logan-GA has strong negative bias. MAI has unacceptably long tails except for the 
case of low noise with long scan duration, i.e. Sc = 1 with 90 minutes scan duration. On the other hand, 
the histogram for Algorithm [1] is close to a Gaussian random distribution; the mean is near zero and the 
distribution is approximately symmetric. Moreover, Algorithm [1] performs well, and is only outperformed 
marginally by MAI for the lower noise and longer time window case. On the other hand, there are some 
situations, particularly for MAI, in which the relative error is less than —100%; in other words, the 
calculated DVs are negative. Such unsuccessful results occur only for the higher noise level (Sc = 2). 
While there was only one such occurrence for the Logan-GA (70 min. with ROI 9) , there were 40 such 
occurrences for MAI, 33 for the shorter time interval of 70 minutes (ROIs 1, 3, 4, 5, 6, 8 and 9) and 7 for 
the longer interval of 90 minutes, (ROIs 1 and 6). The reason for the negative DV for MAI is discussed 
in Section 15.41 From the results for the higher noise Sc = 2 we conclude that Algorithm [T] using the 
shorter 70 minutes scan duration outperforms the other algorithms, even in comparison to their results 
for the longer scan duration. 

Obviously Algorithm [T] is more expensive computationally than Logan-GA and MAI. In the simula- 
tions, the average CPU time, in seconds, per TTAC was 0.00083, 0.00057, 12.2 and 0.0036, for Logan-GA, 
MAI, KA and Algorithm [H respectively. The high cost of the KA results from the requirement to use 
a nonlinear algorithm. Because the KA requires a good initial estimate for the parameters the cost is 
variable for each TTAC; it is dependent on whether the supplied initial value is a good initial estimate. 
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Figure 4: Histograms for normalized error (in percentage), 100(DV C 
four methods. The error ranges are presented in the legends. 
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Figure 5: Histograms for normalized error (in percentage), 100(DV es t — DVt) /DVt, of the results for ROI 3 and four methods. 
The error ranges are presented in the legends. 
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Figure 6: Cs(i)/Cr(i) against time for all test ROIs except ROIs 3, 5 and 11 for the first 90 minutes (a) and 720 minutes 
(b). Dotted vertical lines are plotted at time t* = 35 minutes (a) and t* = 100 minutes (b). The curves for ROIs 3, 5 and 
11 are similar to those for ROIs 1, 4 and 10 resp.. 



Indeed the KA results take from 8 to 25 seconds, while the costs using the other methods are virtually 
TTAC independent. 



5. Discussion 



5.1. Equilibrium Behavior and Dependence on the Size of k^ 

The graphical analysis methods of Logan-type rely on the assumption that the ratio Cs(i) to Cx(i) 
is approximately constant within a chosen window [t* , T] . This ratio is plotted against time for the 



simulated data for ROIs 1 to 11 in Figure 6(a) It is clear that the ratios for ROIs 1, 3 and 6 have not 
reached equilibrium even by 90 minutes. These are the three data sets with the largest bias reported in 
Section [4.21 and with smallest k& (resp. k$). It is certain that equilibrium is eventually reached. These 
curves first increase to a peak at about 120 minutes for ROIs 1 and 3 and at about 180 minutes for 
ROI 6 and then decrease before reaching approximately constant values (Figure [6(b)[ ). On the other 
hand, increasing the scan duration to more than two hours is not practical. Moreover, as illustrated in 
Figure [71 using the linearity of Jq C^{T)dr /Cx(t) versus J* * C^{T)dT /C p {t) to verify whether equilibrium 
has been reached may be misleading. For example, it would appear that all eleven data sets have achieved 
equilibrium after roughly 35 minutes. The arrow in Figure [7J points to the marker corresponding to the 
data calculated at the middle point of the frame from 35 to 40 minutes. 

We illustrate the relation between the bias in the estimate of DV calculated by Logan-GA and k± in 
Figure El As discussed in Section [231 a small value of k^ may cause a large variation in s(i). This graph 
verifies that the magnitude of the bias decreases as k^ increases, further verifying that large bias in DV 
may arise purely due to modeling assumptions in the absence of noise in the data. 



5.2. The effects of quadrature error 

Both Logan-GA and MAI, (|2|) and ((3]) resp., require the calculation of integrals Ct(t)<1t and 

Jo Cp( r )d r - Assume the noise-free measurements Ct(U) are derived from the integral over the ith frame 
duration. Thus we can easily recover its integral without introducing error while quadrature error for 
calculation of Jjj C p (r)dr due to using a limited number of plasma samples is unavoidable. The accuracy 
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Figure 7: f c t C T (T)dr/C T (t) (y-axis) against f* Cp(r)dT/C T (t) ( x-axis) for all test ROIs except ROIs 3, 5 and 11 for the 
first 90 minutes. The last eight points correspond to the time interval 35 to 90 minutes. The curves for ROIs 3, 5 and 11 
are similar to those for ROIs 1, 4 and 10 resp.. The arrow points to the first frame falling in this interval for ROI 6. 
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Figure 8: The bias in the Logan-GA estimation of the DV against the value of k± for the eleven ROIs, assuming noise- 
free data, a scan duration of 90 minutes and t* — 35 minutes. The specific data pairs (&4, bias) are, for ROIs 1 to 11, 
respectively, (0.0106, -0.4231), (0.0230, -0.1024), (0.0106, -0.4318), (0.0170, -0.3286), (0.0150, -0.3874), (0.0050, -1.3459), 
(0.0420,0.0267), (0.0182,-0.1130), (0.0251,-0.0736), (0.0256,-0.0195), and (0.0253,-0.0288). 
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of the numerical quadrature impacts the accuracy of the parameter estimates. Note that we classify the 
noise effects as another source of bias in DV. 

We recalculate the DV for the experiments reported in Section B~2l but now using numerical quadrature 
for calculation of f~ C p (r)dr with data sampled one time point per time frame. The bias for each ROI 
of the estimated DV using 90 minutes scan data with t* = 35 minutes is —11.83%, —2.99%, —11.91%, 
-4.88%, -5.64%, -30.49%, -1.22%, -4.61%, -2.81%, -2.40% and -2.63% when calculated using 
Logan-GA, and -12.02%, -3.10%, -12.10%, -5.01%, -5.77%, -30.42%, -1.28%, -4.87%, -2.93%, 
—2.61% and —2.86% calculated using MAI. It is interesting to note that the DV calculated for ROI 7 is 
no longer an over-estimate. This does not contradict the result of Theorem [TJ which predicts that the DV 
for ROI 7 will be over-estimated due to model error, provided that the other aspects of the calculation 
are accurate. Now using a less accurate quadrature the negative bias due to quadrature error canceled 
the positive bias due to the model error. Indeed, for all eleven test cases the impact of the less accurate 
quadrature is to shift the bias down, i.e. it is more negative as compared to the equivalent more accurate 
calculations shown in Table [2j 

5.3. Bias and classification between AD and NC subjects 

In the eleven simulated ROIs, large under-estimation of the DV calculated by Logan-GA and MAI is 
observed for ROIs 1 (NC Cort), 3 (NC PCG) and 6 (AD Cere). A lower value of the DV in the cortical 
regions of NCs and in the cerebellum for AD subjects will result in under-estimation of the DVR for NCs 
and over-estimation of the DVR for AD subjects when the cerebellum is used as the reference region for 
the DVR calculation. Thus, the difference between AD and NC can be artificially enhanced, and viewed 
as a positive outcome associated with the bias of Logan-GA and MAI. This conclusion, however, can not 
be generalized. It is unknown whether it is always the case that AD/NC have small/large k$ in cerebellar 
regions and relatively large/small k^ in cortical regions. Confirmation of these assertions would suggest, 
based on the discussion in Sections 12.21 and 15.11 that the DVR is over-estimated for AD subjects and 
under-estimated for healthy subjects (also see Figure [8]). In addition, more subtle differences, such as 
the ones between mild cognitive impairment (MCI) and NC, or among NC with differential genetic risk 
for AD, may make the effects of bias much less predictable. Consequently, we evaluate the quantification 
methods based on their bias because the goal of these methods is to estimate the DV as accurately as 
possible. 

5.4. When does MAI fail? 

As noted in Section 14.51 MAI generates some results with negative DVs. Such results are reported 
as unsuccessful in Ichise's original paper [jj. Careful study of these results shows that the negative DVs 
arise when —1/6 has the wrong sign. For most radioligand binding studies 1/6 is a small positive number 
because 6 > (k^ + k&) / {fak^) , which is usually larger than 10, see Remark (fTJ) of Algorithm [TJ Thus 
a small error in the estimate of —1/6 due to large noise in the data may change its sign. This in turn 
impacts the sign of the estimate of the DV. 

6. Conclusions 

In this article, we quantified the model error in estimating distribution volume using graphical analysis 
methods. We described the conditions under which the DV is either over- or under-estimated, and 
quantified the bias caused by model error. We validated our findings through simulations with noise-free 
data. To reduce the impact of model error, we added a simple nonlinear term to the fundamental linear 
model MAO, and presented a new algorithm for its solution. Simulations with noisy data demonstrate that 
the new algorithm is cost-effective and robust even for shorter scan durations. For PIB-PET studies, the 
new method using shorter scan data (70 minutes) outperforms, or is at least as good as, Logan-GA, MAI 
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and KA methods using longer scan data (90 minutes). The proposed approach can be easily extended 
for DVR estimation. This is a focus of our future work. 
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8. Appendix A: fundamental theory for Corollary Q] 

Here we present the theoretical result from which Theorem [T] is obtained. We use the notation that 
a = (ai , Gt2, ■ • • , a n ) T and b = (b\, 62, • • • , b n ) T , are vectors with entries a% and hi, resp. The notation 
a/b and aob denotes component wise division and multiplication, namely entries ai/bi and a^, ||a||i is 



Y^i=i l a *l an d Il a ll2 = yja% + a| + ■ ■ ■ + a£ is the Euclidean norm. We call a decreasing (increasing) if 
a i > 02 > • • • > On ( a i < a 2 < • • • < a n)> an d non-constant decreasing (non-constant increasing) if it 
is decreasing (increasing) and at least one of the > (<) signs is strict, > (<). If all of the > (<) signs 
are strict, we call a strictly decreasing (strictly increasing). A vector a is constant if = a for some 
constant a and for all i. 

Lemma 1. ( Chebyshev's sum inequality [2^] ) Given real numbers 
0.1 > d2 > • • • > a n and bi > 62 > ■ ■ ■ > b n , then 



it^y^taMltbA. (22) 
k=l \ k=l / \ fc=l / 

Similarly, if a\ > 02 > ■ ■ ■ > a n and b± < 62 < • • • < b n , then 

;E»As(^"')(^4 (23) 

k=l \ k=l J \ k=l / 

In the above Chebyshev's sum inequalities the numbers are not required to be positive and the equality 
is true if and only if one of the two vectors, a or b, is a constant vector. If a and b are positive vectors, 
the Chebyshev's sum inequalities can be expressed as a r b > 1 1 a| 1 1 1 1 lo 1 1 1_ and a T b < - ||a||i||b||i. 

Lemma 2. Ifp, q and s are positive real vectors, of which p is a increasing vector and q is a decreasing 
vector, then 

1. ||q|||p s — p T qq T s > i/s/q is a non-constant increasing vector. The inequality is strict ifp is 
strictly increasing. 

2. ||q||2P T s — p^qq^s < i/s/q is a non-constant decreasing vector. The inequality is strict ifp is 
strictly increasing. 

3. ||q|||p s — p T qq r s = i/s/q is a constant vector, 

4. ||p|||q T s — p T qp r s > if s/p is a non-constant decreasing vector. The inequality is strict ifp is 
strictly increasing. 

5. HpH^q^s — p T qp r s < if s/p is a non-constant increasing vector. The inequality is strict ifp is 
strictly increasing. 
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6- ||p|| q s — p qp s = i/s/p is a constant vector. 

Proof. We only prove the first case. The proof for the other items follows similarly. We use mathematical 
induction. For the lowest dimension n = 2, 

l|q|| 2 p s-p qq s 

= (Qi + ?i)(Pl s l +P2S2) - (piqi +P2Q2)(qiSi + Q2S2) 

= qlP2S2 + qlp\S\ - Piqiq2S2 - £>2<?2<?l.Sl 

= (91 S2 - g 2 si)(gip 2 - 

= (91P2 - 92Pi)((gig2)(— - — )) 

?2 gi 

> 0. 

The last reduction follows from the monotonicity of p, q, which implies q±p2 — q2P\ > 0, and the non- 
constant increasing assumption of s/q, which guarantees ^ — ^ > 0. When p is strictly increasing 
Q1P2 — q2Pi > 0. Under this condition ||q||2p T s — p T qq T s > for n = 2. Assuming the inequality 
||q||2P r s — p T qq T s > is true for dimension n = i, i.e. 

i i i i 

^2 Qk^PkSk - ^puqu Y^j ikSk > 0, 

k=l k=l k=l k=l 

then for n = i + 1 

11 \i2 T T T 

NI2P s ~p qq s 



= (%2<it + qi+i)C^2pk s k + Pi+iSi+i) - C^Pkqk + Pi+iqi+i)(^2qkSk + qi+is i+1 ) 

k=l k=l k=l k=l 

i i it 

= C^ql^PkSk -^Pkqk^qkSk) 

k=t k=l k=l k=l 

iii i 

+{p i+1 s i+ i ^qj- qi+iSi+i^Pkqk) + (qf+i ^JpfeSjfe - Pi+m+i ^2 qkSk ^ 

k=l k=l k=l k=l 

i i 

> + Sj+i ^ qk{qkPi+i - Pkqi+i) + ^2 s k (q i+1 p k - p i+ iq k ) 

k=l k=l 

i 

= ^iqkPi+i - Pkqi+i)(qkSi+i - qi+is k ) 



k=l 

i 

Si+1 S k . 



= J2{(qkPi+i-Pkqi+i)(qkqi+i)) 





k=l qi+1 qk 



The last reduction is based on the monotonicity of p, q and s/q. When p is strictly increasing q k Pi+i — 
PkQi+i > for all k < i the inequality will be strict because at least one of the terms — ^,k = 1, • • • ,i, 
is positive based on the monotonicity condition. The result thus follows by induction for all integers 
n > 2. □ 

The following corollary now follows immediately by observing that s/q increases when s increases and 
s/p decreases when s decreases. 
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Corollary 2. Ifp, q and s are positive real vectors, of which p is a strictly increasing vector and q is a 
decreasing vector, then 

1. ||p||2q T s — p T qp r s > if s is a decreasing vector. 

2. ||q|||p s — p T qq T s > if s is an increasing vector. 

Lemma 3. If p, q, r and s are positive real vectors, of which p is strictly increasing, q is decreasing, 
and p, r, s and x* satisfy px* — s = r; and [x, b] = argmin||px — 6q — r 1 1 1 ; then 

1. £/ie estimated solution x and exact solution x* are related by 

• x > x* i/s/q is a non-constant decreasing vector, 

• x < x* i/s/q is a non-constant increasing vector, 

• x = x* if s/q is a constant vector; 

2. the following inequality is true without any monotonicity assumptions: 

|f - ,|S l|p|IIINII-t^ F(g) - (24) 

3. the sign of the intercept b is determined as follows: 

• b > if s/p is a non-constant decreasing vector, 

• b < if s/p is a non-constant increasing vector, 

• 6 = if s/p is a constant vector; 

4. given x = x* , the LS solution ofpx — 6q ~ r /or b is b = q T s/||q||2,' 

5. given b = q T s/||q|||, the LS solution ofpx — 6q ~ r for x and the true solution x* have the same 
relationship as stated in the first conclusion of this theorem. 



Proof. It is easy to verify that the LS solution of px — 6q ~ r is 

ii 1 1 9 T 1 7 1 T 1 ii 1 1 9 T 1 T 1 T" 

q oD r — p qq r - — P 2Q r + P QP r 
x = — — — b = — - — — 

llpllillqlll - (p T q) 2 ' llplllllqlli - (p T q) 2 

The proof then follows as outlined below: 

1. Replace r in the expression for x, with px* — s. Then 

l|q|l2P T (P x * ~ s ) ~ p T qq T (px* - s) 

llpllillqlll - (p T q) 2 

T T II 112 T 

* p qq s — IIQH2P s 

I|p|l2lq|l2 - (p q) 

and the results immediately follow from Lemma[2] dU)-© and the fact ||p||2||q||2 > P T q when p is 
not linear proportional to q. 

2. Because 

T T 11 n2 T 
p qq s - ||q|| 2 p s 

T / \T — 11 n2 / \T — 

= p q(q°q) s- l|q|| 2 (p°q) s 

< P^qllqlli max(s"i) — ||q||2P T q ' min(si) 

i i 

= P T q||q|l2( max (^) - min(si)), 

% 1 
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and similarly 

p T qq T s - ||q||2p T s > p T q||q||2(min(si) - max(sj)), 

We have 

|p T qq T s — Hqll^p^sl < p" r q||q||2(min(sj) — max(sj)). 

i i 

Using the fact HpIIHIqIII ~~ (p T q) 2 > an d p5|) . we conclude the inequality is true. 
3. Again we replace r with px* — s, then the expression for b becomes 

- — ||p||2q T (px* - s) + p T qp T (px* — s) 

llpllllqlll - (p T q) 2 



i 1 1 9 T T T 

Ip^q s-p qp s 
llplliNI! - (p T q) 2 ' 



(26) 



The results immediately follow from Lemma[2] flU)-© and the fact Hplbllqlb > P T q when p and q 
do not have the same direction. 

4. This result is easily verified. 

5. Given b = q T s/||q||2, the LS solution of px — 6q « r for x is 

£ = n 1 i9 P T (q^ + r ) 
WpM 

i i t q Ts , t i * \\ 
= 7i — (p qTnT2+P (p x ~ s )) 
IIpIIs Ilqll2 

T T ii ii 9 T 

* P qq's- NIIp^s 

' ii 1 1 9 II 1 1 9 

The results now follow from Lemma [2l 

□ 

We now transform the exact equation to p/qx* — s/q = r/q and rewrite the results using vectors 
P = p/q> s = s/q and f = r/q. Correspondingly, we find the LS solution of px — eb ~ f for e = (1, 
1,-,1) T - 

Corollary 3. Ifp, f and s are positive, of which p is strictly increasing, p, f, s and x* satisfy px*— s = r; 
and [x, b] = argmin||px — be — f |||, then 

1. the estimated solution x and the exact solution x* are related by 

• x > x* if s is a non- constant decreasing vector, 

• x < x* if s is a non- constant increasing vector, 

• x = x* if s is a constant vector; 

Moreover, the following inequality is true without any monotonicity assumptions. 

\x-x*\ < „ „ 9 y(s). (27) 

n||p|l2 ~ IIpIIi 

2. The sign of the intercept b is determined as follows: 

• b > if s/p is a non-constant decreasing vector, 

• b < if s/p is a non-constant increasing vector, 
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• 6 = i/s/p is a constant vector. 
In addition, 

n 

• b > ^^Sj/n if s is a non-constant decreasing vector, 

i=i 

n 

• b < ^^Sj/n if s is a non-constant increasing vector, 

i=l 
n 

• b = Si/n ifs is a constant vector; 

i=i 

n 

3. Given x = x* , the LS solution of px — be ~ r for b is b = s-Jn; 



i=l 



4. Given b = ^^Si/n, the LS solution of px — be m r for x and the true solution x* are related 



as 



i=l 

stated in the first conclusion of this theorem. 

Proof. Most results are a direct Corollary of Lemma [3] by setting q = e. We only prove the new results 
P and ©. 



1. We just need to prove the bounds for \x — x*\. Setting q = e in (|25p we have 

IIpIIi ~ n P T s 

x = x* + — i= „_,.- . (28) 

"IIpII^ - IIpIIi 

Because 

IIpIIi / * 8% ~ np T s < n ■ max(s.j) ||p||i — n ■ min(sj) ||p||i = n||p||i(max(sj) — min(sj)), 

i 

IIpIIi / ^ &i — np T s > n ■ min(sj) ||p||i — n ■ max(sj) ||p||i = n||p||i(min(sj) — max(sj)), 

^ — J i i i i 

i 

and nllpUl — ||p||i > we obtain 

I ||p||i y^Ji - np T s\ 



\x — x 



HIpIII - IIpIIi 

n||p||i(maxj(5j) - min^)) 

HIpIII - IIpIIi 
"IIpIIi T/ r/ s x 
"llplla - IIpIIi 

2. Setting q = e in (j26|) we have 

t llPllil|s||i - ||p||ip t s 

b = [73772 \TU2 ' 

»P - Pi 
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The results on the sign follow from Lemma [3] ([3]). For the remaining three inequalities, we only 
prove the case for which s is decreasing. Proofs of the other two are similar. Setting q = e in f|26[) 
we have 

t ||p|||||s||i - ||p||ip T s 



> 



n \\P\\l - IIpII 2 
IpjlflNi - IIpIIi^IIpIIiI 
"IIpIII - IIpII? 

n 

»||1 _ i=] 
n n 



□ 



9. Appendix B: component-wise perturbation analysis for LS solution of (|18|) 

In Remark 2, we claimed that "the estimate of DV is much more robust to noise in the formulation 
than are the estimates of A and B because Jq C p (r)dr is much larger than both Cx(i) and Cs(i) for 
t > t* n . Here we present a theoretical explanation, which is helpful for algorithm design in quantification. 
Instead of considering a general linear equation, which is out of the range of this paper, we assume a 
system of equations Ax. = y with only two independent variables x = [x\, X2[ r ■ The two columns of the 
system matrix A are denoted by ai and a2, i.e. A = [ai,a2]. 

Theorem 1. Suppose the linear system Ax. ~ y + e, for A = [ai, a.2\, has the exact solution x = [x^x^], 
the uncorrelated noise vector e obeys a multi-variable Gaussian distribution with zero means and common 
variance a 2 and that \\sli\\ » \\&2\\- Then least squares solution x = [xi,X2] T has the following statistical 
properties 

1. E(xi) = x\ and E{x2) = x\, and 

2. Var(xi) << Var(x 2 ). 

Proof. We assume matrix A has the following singular value decomposition 

( cos6 s ' m0 \ (00) 
\ -sm9 cos9 J ' K ' 

in which s\ > S2- Then 

x = VS^U T (y + e) = x* + VS*U T e, 

where 

ot_ / 1/si ... 0\ 

V l/8 2 ••• J' 

Because U is an unitary matrix and ||ai|| >> ||a2|| we immediately derive the the following inequality 
from equation (|29|) : 

si cos 2 e + si sin 2 e » s\ sin 2 + s 2 , cos 2 9. 



A= [ai,a 2 ] = USV 1 



U 



( si 






S2 



\ 
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LS solution histogram: x 1 and x 2 , ( 1% noise in data.) 
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This inequality is equivalent to (s 2 — s 2 ,) cos 2 9 + s 2 , » (s 2 — s 2 ) sin 2 9 + s 2 , which implies cos 2 9 » sin 2 0, 
i.e. cos 2 9 m 1 and sin 2 flwO, and sf >> s 2 . If we denote the two rows of matrix VS^ by qi and q2 than 



|qi|| 2 = sin 2 9/s\ + cos 2 9 j s\ 
|q 2 || 2 = cos 2 9/s\ + sin 2 9 j s\ 



1/s 2 + cos 2 6>(l/s 2 , - l/sf), 
l/s\ + sin 2 9{l/s 2 2 - l/s\). 



Because cos 2 9 » sin 2 9 and l/s| >> 1/s 2 we conclude ||q2|| 2 >> ||qi|| 2 - If we let pi and p2 be the two 
rows of matrix VS^U T then ||pi|| 
Let 



|qi|| and ||p2|| = ||q2|| because U is unitary. Thus ||p2|| IIPil 
d = x-x* = VS*U T e. 



It is clear E[d\) = and E{di) = because the means of e are zero, and Var(di) = YliPii a2 = HPilT " 
and Var(d2) = YliPli' 7 ' 2 = IIP2|| 2(j2 resp.. Therefore Var(di) << Var(d2). Because d = x — x* we 
conclude E(x) = x* and Var(x\) « Var(x2) □ 

This result is illustrated by the following simple example: 




The first column is much larger than the second column. If we add 1% noise to the right hand side, i.e. 
e\ ~ N(0, 0.05), £2 ~ iV(0, 0.09) and £3 ~ N(0, 0.115), and perform simulation with 1000 realizations the 
distribution of the resulted x\ and X2 are illustrated in Figure [9l These results are consistent with the 
conclusions in Theorem [TJ 
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10. Appendix C: derivation for equation (1121) 

Integrating © and ([6]) from to t we obtain 

C F+NS (t) = Kx f C p (T)dT-(k 2 + h) I C F+NS (r)dr + k 4 f 
Jo Jo Jo 

C S (t) = k 3 I C F+NS (T)dT - h [ C s (r)dr, 
J o Jo 

= h C F +Ns(T)dT - k± I (C t (t) - C F+NS (T))dT, 

Jo Jo 

= -k 4 / C T (r)dT + (k 3 + k A ) / C F+NS {r)dr. 
Jo Jo 

Taking the sum of equations ([30]) and ([3T]) yields: 

Cr(t) = K x f ' C p {r)dr - k 2 f C F+NS (r)dr, 
Jo Jo 

and canceling J * C F+ Ns(i~)dT from ([32]) using ([33]) gives: 

C s (t) = -h J* C T (r)dr + jT* C p (r)dr - C T (t) 

This can be transformed to f) 121) immediately by using = ^(1 + f^)- 
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